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I. INTRODUCTION 
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(Date: February 1, 2008) 

We consider the Brownian motion of a quantum mechanical particle in a one-dimensional 
parabolic potential with periodically modulated curvature under the influence of a thermal heat 
bath. Analytic expressions for the time-dependent position and momentum variances are compared 
with results of an iterative algorithm, the so-called quasiadiabatic propagator path integral algo- 
rithm (QUAPI). We obtain good agreement over an extended range of parameters for this spatially 
continuous quantum system. These findings indicate the reliability of the algorithm also in cases 
• for which analytic results may not be available a priori. 
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Exactly solvable systems have a special status among physical models. Although oversimplified in many cases, they 
. may serve as starting point for testing the reliability of methods which can then be transferred to more realistic, but 
only numerically solvable models. An important class of such models are quantum systems coupled to a dissipative 
, environment and being driven by a time-dependent external field . A wide variety of physical phenomena have been 
■ described by this kind of models, e.g. electron Q| and proton [|j transfer, tunneling processes of a macroscopic spin 
, hydrogen tunneling in condensed phases |^ , single defect tunneling in mesoscopic quantum wires , or tunneling 
"Y^ of the magnetic flux in a SQUID j^], to name but a few. Usually such dissipative quantum systems consist of a model 
Hamiltonian bilinearly coupled to a bath of harmonic oscillators. Additional external time dependent driving fields 
^ render the mathematical solution even more difficult or even impossible. 

One of the few analytically tractable time-dependent dissipative quantum systems is the parametrically driven 
harmonic oscillator whose analytic solution was found by Zerbe and Hanggi in Ref. |^ . A physical realization of this 
model is the Paul trap j^], which provides an oscillating quadrupole potential for the enclosed ion. Furthermore, the 
^ ] parametrically driven dissipative harmonic oscillator may serve as a benchmark for approximation schemes which were 
developed for more general dissipative systems . The interesting feature is that the parametric driving induces a 
non-trivial quasienergy spectrum [Q, in contrast to additive driving where the quasienergy spectrum coincides with 
the spectrum of the undriven system apart from a constant shift. This is further corroborated by the fact that the 
solution of the parametrically driven linear oscillator can be utilized to obtain solutions of certain nonlinear dynamical 
systems [ pd| . 

Powerful approximative numerical procedures for simulating dissipative and possibly time-dependent quantum sys- 
tems are the Quantum Monte Carlo method |l2| and the quasiadiabatic propagator path integral algorithm (QUAPI) 
developed by Makri and Makarov |p^ . The former method works very well for problems involving path integrals 
in imaginary time, however, the calculation of real time path integrals is afflicted by the so called sign-problem due 
to the rapidly oscillating integrand. The QUAPI algorithm has been applied to low dimensional dissipative systems 
such as the driven spin-l/2-particle coupled to a harmonic oscillator bath (driven spin-Boson-system) p3| and to 
the driven double- well potential in order to study quantum hysteresis and quantum stochastic resonance [ ]l^ , p^ . 
Moreover, the QUAPI algorithm has recently been used as a basis for a very efficient memory equation algorithm for 
spin-boson- models ]l6| |. 

The purpose of this paper is to apply the QUAPI algorithm to the parametrically driven harmonic oscillator and to 
compare the results with the analytic solution from Ref. |^ . Whilst harmonic oscillator systems are known to exhibit 
some untypical features this is not the case with respect to the QUAPI algorithm. Our results thus show that not 
only intrinsically discrete models like the spin-boson-system but also spatially continuous systems can be accurately 
described by few energy eigenstates if the temperature is restricted to a moderate regime. Most importantly, this 
is the first work in which the numerical approximative QUAPI results are compared against analytic solutions of a 
spatially continuous driven dissipative quantum system. 
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The paper is organized as follows: In section ||, we introduce our model of the parametrically driven dissipative 
quantum harmonic oscillator and briefly review the analytic solution given in Ref. 1^. Section III is devoted to a 
short review of the QUAPI method. The comparative main results are presented in section IV, before we give the 
conclusions in section 



II. THE MODEL AND ITS ANALYTIC SOLUTION 



In this section we briefly review the analytic solution of the parametrically driven dissipative harmonic oscillator 
from §. 

A quantum particle with mass M , position operator x and momentum operator p moving in a one-dimensional 
harmonic potential with periodically modulated curvature is described by the Hamiltonian 

2 

(1) 

Following the common approach [Q to model the influence of the environment by an ensemble of harmonic oscillators, 
the bath Hamiltonian Hb (including the interaction with the system) is given by 



He ^ E H, (x, = 5: i[£L + __iL^3c)]. (2) 



] J J J J 



and the whole system is described by the Hamiltonian H(i) = Hs(t) + Hb- In the case of a thermal equilibrium bath, 
it turns out that its influence on the system is fully characterized by the spectral density 



JM = |E;|;-^(--.)- (3) 

With the number of harmonic oscillators going to infinity, we arrive at a continuous spectral density. In the following, 
we choose for the sake of definiteness a truncated Ohmic spectral density, i.e. 

J{uo)^M-tuf,{uo,uo,). (4) 

Here, 7 is the coupling strength to the heat bath and fc{Lu,LUc) denotes a cut-off function which avoids unphysical 
divergences due to high-frequency bath modes. For our calculations, we consider two examples for the cut-off function: 
(i) a smooth exponential cut-off 

fc{uj,uJc) = e^p{-uj/ujc) (5) 

and (ii) a step-function 

/c(w, Wc) = 6(wc - t^) (6) 

with cut-off frequency tOc ^ too,^ (see discussion given below). 

We choose a factorizing initial condition of Feynman- Vernon form jl^ which means that at time t = to, the full 
density operator W(io) is given as a product of the initial system density operator Pg(to) and the canonical bath 
density operator at temperature T = l/fcs/S, i.e., 

W(<o) = Ps(^o) ^B ' exp(-/3HS^) , (7) 

where Z^^ = Tr exp(-/3H|^) and 

By way of integrating out the bath degrees of freedom in Eq. (||) one obtains the following one-dimensional Heisenberg 
equation for the position operator x, i.e. 

/•* 1 
x(i) + / l{t- t')±{t')dt' + (co2 + e cos nt)^{t) = — r(t) - 7(i - to) x(io) , (9) 
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with the friction kernel given by 



2 /■°° , J{uj) , , 



r(t) is a time-dependent fluctuating (operator) force 
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Pj(^o) 



s\n{uij[t - to)) + qj{to) cos{ujj{t - t^)) 



(10) 



(11) 



which contains the initial conditions of the bath and of the particle's position at time to. The last term on the r.h.s. 
(proportional to x(io)) in Eq. (||) is the so-called initial slip, caused by the specific choice (^ of the initial conditions. 

Exploiting the thermal distribution of the bath one recovers the usual connection (via J{uj)) between the random 
and the frictional forces of the bath in Eq. (O) in the form of the fluctuation-dissipation-relation, reading, t > t', 



(r(t)r(i'))/3 = Tr [Zb 1 exp(-/3H° ) r(i)r(i')] = hL{t - t') 



(12) 



1 f°° 

L{t) = - I dujj{uj) 
^ Jo 



coth ^ cos(ijjt) — ism{ujt) 



(13) 



where the subscript /3 indicates thermal averaging performed with the canonical density operator for Hg defined in 
Eq. (^. The response function L{t) will play an important role in the numerical QUAPI algorithm. 

It turns out |^ that for the description of the parametric dissipative quantum oscillator the solution of the classical 
deterministic limit {h ~> 0,T 0) with uj c ^ oo plays a prominent role. Thus, in Eq. (||) the position operator x 
is replaced by the classical coordinate x and J^^ ^{t — t')x{t')dt' goes over into jilt). Moreover, on the right hand 
side of Eq. (|9|), the fiuctuations r{t) are zero and the initial slip is also omitted, which can be achieved by either a 
somewhat different choice of the initial conditions than in Eq. (7|) or by replacing the coupling coefficients cj in Eq. 
(^ by Cj Q{t — t^) so that Hb and Hg from Eq. (^) coincide at t = to. For convenience wc furthermore introduce 
scaled quantities 



i = §t, i(i) 



y/Am/2hx{t = If), Cjo = ^LUO, 



^-7 



T 



2k_B_rp 

nn ^ ' 



UJc 



(14) 



In the remainder of this paper, we exclusively use dimensionless quantities but omit all the tildes for the sake of better 
readability. In order to recover the dimensionful quantities, one has to re-introduce tildes wherever it makes sense 
and then exploit Eq. (1^). By substituting x{t) = y{t) exp[— 7(t — to)/2] we arrive at an undamped oscillator equation 
for y which is the well-known Mathieu equation 



y(t) + (cj^ - i- + 2ecos2t)y{t) = 



(15) 



Its mathematical properties like stability and instability regions in the parameter space are well known p8| . Never- 
theless, there exists no closed analytic expression for the solution and the equation has to be integrated numerically. 
In the following, we will need two linear independent solutions ^i(t),i = 1,2, of Eq. ( p^ belonging to two different 
sets of initial conditions 

$i(to) = 0, $i(to) = 1, 

$2(to) = 1, $2(to) = 0. 

They can be determined numerically, e.g. by means of a regular fourth-order Runge-Kutta integration of the Mathieu 
equation (p^. 

Let us return to the dissipative quantum parametric oscillator. The quantities of interest are the variances of the 
position and the momentum operator, i.e. 



a.4t) ^ (x^(t)) - (x(t))^ 
a^pit) = i(x(t)p(t) +p(<)x(t)) - (x(t))(p(t)), 
a,,it)^{p'{t))-{pit)r. 



(16) 
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Here, the quantum mechanical expectation value is understood as usual as (•) = Tr[p{t)-]. By determining the 
propagator U(i,to) = Texp(— ? j^^ dt'H{t')/h) (T is the time ordering operator) for the driven dissipative system 

according to the reduced density matrix p{t) = rrBQt;i(U(i, fo)W(to)U~-^(t, ^o)) can be calculated analytically. 
Here, W(io) denotes the full density operator at time and TrBath the trace over the bath degrees of freedom. 
Having obtained the reduced density operator p(i), the quantum mechanical expectation values in Eq. [T^ can be 
evaluated. After some algebra, we find for the dimensionless variances the expressions 



e-T(*-*«){[*2(i) 



7, 



■^i{t)fal^ + 2$i(t)[$2 W - J$iW]cT°p + <^l{t)al^] + S,,(t), 



(t) + [cj^ + 2ecos(2i)](T,,(t) - Epp(i). 



(17a) 

(17b) 
(17c) 



Thereby, we have rectified p9[ | some minor misprints in |^ and simplified the equations in |8| for a^p and cTpp. Here, 



denote the initial variances of the uncoupled system at time t — which depend on the choice of 
the initial state for the bare system Hs(to)- The initial conditions for Eqs. ( |l7| ) at time t = ^o" are given by 



' XX ' '^xp S'Hd (Jpp 



l^lx 



xp 1 



<^pp{tt)^l^<x-'^l<p + ' 



pp- 



(18) 



The discontinuity of the variances at time is a well-known consequence of the initial slip term in Eq. (^; it is 
due to the factorizing initial condition (|^) . The first terms in the three equations ( p7| ) posses the same form as in the 
classical case. The specific quantum mechanical features enter via the functions Yixx{t) and Spp(<), which read 



7 / ^ \ \ I 7 

T.x,At) = -j dwtj/c(w,Wc) coth(— )< / dsG(t, s)exp(-(t - s))cos(ws) 



ds G{t, s) exp( — (t — s)) sin(cjs) 



T,pp{t) - dujujfc{uj,u;c)coth{—) ds G(t, s) exp(-(i - s)) cos(w(t - s)) , 



'2T 



(19a) 
(19b) 



to 



where G{t,s) = $i(i)<I>2(s) — $i(s)$2(i)- While in Eq. (|l9| ) a general form of the cut-off function fc{i^,ujc) with 
ujc 3> luq is kept, the analytic solution ( p^ is based |^ on the assumption of a strictly Ohmic classical dynamics 
{ujc —> oo) in Eq. (|l0|). The consequence of this assumption is the discontinuity at t ~ to in Eq. ( p^ when the 
system-bath-interaction is switched on instantaneously. A finite cut-off in the spectral density J{uj) in the damping 
kernel (|l^) would induce a smoothend time evolution of the variances ([l7| ) close to i = io on a time scale w^T^. 

The relations ( p^ , pj| ) are evaluated by standard numerical methods. The efficiency is improved if one applies 
Floquet's theorem for the fundamental solutions ^j{s). Then, the periodic part of the Floquet solutions can be ex- 
panded in a Fourier series and the integrations over the intermediate times s in Eq. (|l9|) can be performed analytically. 
Finally, the remaining w-intcgrations and the sum over the Fourier modes can be readily carried out. 



III. NUMERICAL SOLUTION WITH REAL-TIME PATH INTEGRALS 



In the following section, we recapitulate the essentials of the QUAPI algorithm. Further details can be found in the 
original works by Makri and Makarov . In order to describe the dynamics of the system of interest it is sufficient 
to consider the time evolution of the elements of the reduced density matrix which reads in position representation 



U{tf,to) = Texp <{ ~i/h / n{t')dt' 

J to 



(20) 



Here, T denotes the chronological operator, W(to) the full density operator at the initial time tg and TrBath the 
partial trace over the harmonic bath oscillators qj. Due to our assumption that the bath is initially at thermal 
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equilibrium and decoupled from the system, W(<o) becomes the product of the initial system density operator Pg(to) 
and the canonical bath density operator at temperature T, see Eq. (^. Then, the partial trace over the bath can be 
performed and the reduced density operator be rewritten according to Feynman and Vernon as 



p{xf,x'f,tf) = J dxodx'^G{{xf,x'j:,tf;xQ,XQ,to)p{xo,x'Q,tQ), (21) 

with the propagator G given by 

G{xf,x'j:,tf;xo,XQ,to) = J VxVx exp {Ss[x] - Ss[x'])^ Tfv[x,x']. (22) 

Ss[x] is the classical action functional of the system- variable x along a path x{t) and J-fv[x, x'] denotes the Feynman- 
Vernon influence functional 

TFv[x,x']^expl^-^ J^^ dt J^^ dt' [x{t') - x'{t')][T]{t' - t)x{t) -r]*{t' -t)x'{t)]j , (23) 



with the integral kernel 



T]{t) ^ L{t) + iS{t) - / —J{lj) (24) 



and the autocorrelation function L{t) given in Eq. (|l3|). As usual, the restriction to paths that satisfy the boundary 
conditions xo{to) = Xq, Xf{tf) = Xf and similarly for x'(t) is understood implicitly in Eq. (|2^). Likewise, the 
dependence of the density operator p in Eq. (21) on the initial time to and on Ps(^o) has been dropped. 

To make the equations numerically tractable, we discretize tf — into N steps At, such that tfc = + and 
split the full propagator over one time step \][tk+i,tk) in Eq. ( pO| ) according to the Trotter formula symmetrically 
into a system and an environmental part: 

U(ffe+i,ife) «exp(-iHBAV2ft)Us(tfe+i,tfc)exp(-iHBAV2ft) , 



Us(ifc+i,ife) =rexp|--y dt'Usit')^ . (25) 

The symmetric splitting of the propagator in Eq. ( p5| ) causes an error proportional to At^ . This error will be studied 
in detail in section ^ below. The short-time propagator Us of the bare system is numerically evaluated by means 
of a Runge-Kutta scheme with adaptive step-size control. Exploiting the approximation (p5|), the propagator in the 
position representation now factorizes as 

{xIV,q,\V{tk+iM)Wli,q',) « {x\[Js(tu+i,tu)\x') J|(5^.|e-^H,(.)At/2/.g-^H,(.')At/2fi|^,^ ^ (26) 

where the Hj(x) are defined in Eq. (H). By exploiting this approximation and performing the partial trace over the 
bath modes in Eq. (pO|), one recovers Eq. (|2^), but now with a discretized version of the propagating function (|2^), 
i.e.. 



p{xf,x'j:;tf) — J dxn ■ ■ ■ J dxN j dx'f^ . . . j dx'pf S{x'j — x'j^) S{x j — xjy) 

x{xN\Us{tf,tf - At)\xN-i) . . . (a;i|Us(io + At,to)\xo) 
x{xo\Ps(.to)K){x'o\lJs\to + At,to)\x[) . . . {x'^^,\U^\tf,tf - At)\x'^) 



xJ^I^{xo,x'q,... ,xn,x'n) . (27) 

Here, J-^py {xq, x'j^) is the discrete Feynman- Vernon influence functional ( p3|) where the paths x{t) and x'{t) consist 
of constant segments Xk and xj., respectively, within each time interval tk — ^At < tk < tk + ■^At and can be rewritten 
in the form 



[_ k=Ok'=k 



Xk' - x',^,][i]k'kXk - Vk'k^'k] \ ■ (28) 



5 



The coefficients {rjk'k} are closely related to their continuous time counterpart ri{t) in Eq. (|2J). Their explicit form is 
lengthy and not very illuminating for our purposes; their detailed form can be looked up in Ref. fl^ . 

To make further progress, it is necessary to approximately bre ak the influence kernel !F'^{xo, in Eq. (|2^ ) 



into smaller pieces. To this end, Makri and Makarov use the fact [gO 13 1 that the real part of the integral kernel Lfi{t) 
typically exhibits a pronounced peak at t = 0, and quickly approaches for t ±oo. The decay to zero depends 
naturally on the choice of the cut-off function fc{uj,ujc), see Eq. (^). This suggests the truncation of r]{t) after a 
certain number K of time steps At and, correspondingly, to neglect rjk'k ii k' > k + K, i.e. 

N inin{N,k+K} . ^ . 

T^Jl\xQ, ...,x'n) « Y[ "Ji^^f"' ~ ^'k']iVk'kXk - Vk'kX'k] > ■ (29) 

fc=0 k'=k ^ 

In doing so, we approximate L{t) by zero for t > KAt, cf. Eqs. ( p^ , p^ . Of course, this truncation induces an 
error in the final result which has to be handled with care. The error becomes increasingly less important for 
increasing temperatures since then, the bath-induced correlations fall off increasingly faster. In other words, for 
higher temperatures the width of the response function L{t) decreases. In the other limit of decreasing temperature 
however, the number K of relevant time-intervals is increasing and in the limit of zero temperature T = 0, it is well 
known Q| that the response function L{t) falls off only algebraically for t — > ±c». Nevertheless, we will see that this 
approach allows to deal with quite low temperatures and produces qualitative agreement with analytic solutions. 

The next goal is to approximate the spatially continuous integrals in Eq. ( p7| ) in terms of finite sums. To this end, 
Makri and Makarov perform a transformation into a basis given by the energy eigenstates \(j)m) of the bare system 
Hamiltonian Hs(ir) (0), but with the driving term clamped to an appropriate but fixed (!) reference time tr, i.e. 

HsiUMm) ^ E„,\(j)„,) , m-1,2,... . (30) 

Em denotes the energy eigenvalues of the static system Hamiltonian Hs(ir)- In certain cases, symmetry properties 
suggest the choice of an appropriate tr- Here, we choose the unperturbed harmonic oscillator as a reference configu- 
ration. This means for our choice of driving to use tr — 7r/4, so that cos(2ir) = in Eq. (|l5|). Reintroducing now the 
thermal bath but restricting ourselves to small-to-moderate temperatures T, the thermal occupation of high energy 
levels Em is expected to be negligible. This argument suggests that the \(j)m) provide a well adapted basis admitting 
a fast convergent truncation scheme. In other words, we may approximately project the dynamics onto the Hilbert 
subspace spanned by the first few energy eigenstates \4>m)i m = 1, M, corresponding to an approximate decomposi- 
tion of the identity operator I « Before doing so, we perform one more unitary transformation within 
that M-dimensional Hilbert space such that the position operator becomes diagonal [discrete variable representation, 
DVR Q]: 

M 

\Um) = ^ Rmm'\(l}m') , {Um\x.\Um') = Xm^'^Sjnm' , m, to' = 1 , . . . , M . (31) 
m' — l 

Exploiting the approximate decomposition of the identity I « 'E,m=i\um) {um\ and the truncation of the bath-induced 
correlations in Eq. (p9[), it is a matter of straightforward but tedious manipulations - starting from Eq. (p7| ) - to 
arrive at the final form of the QUAPI recursion scheme. In particular, the integrals in Eq. (|2^) turn into finite sums 
due to the transformation ( ^l|) into the DVR-basis. We do not present the detailed form here and refer the reader 



again to the original literature |13 



The above introduced restriction to a finite dimensional subspace induces an error in the evaluation of the reduced 
density matrix. However, as we will also discuss below, this error behaves in a controlled way if the relevant parameters 
such as the temperature and the damping are chosen in a moderate regime. This means that for increasing temperature 
increasingly more DVR-states are necessary to describe the dynamics appropriately. Note that in this regime however, 
the number K of the relevant memory time steps is decreasing. In the opposite limit of decreasing temperature, the 
number M of relevant basis states can be chosen rather small. In this low-temperature limit the number K of memory 
time steps can therefore be increased. Moreover, we note that the restriction of the dynamics (at long times) to the 
M-dimensional Hilbert subspace is not allowed for systems with an inherent diverging dynamics. This is also seen 
in our example of the parametrically driven dissipative quantum harmonic oscillator for a parameter choice in an 



instability region of the Mathieu equation (n5[) , see section IV below 



The efficiency of the QUAPI algorithm is based on the choice of the two free parameters M (the number of basis 
states) and K (the length of the memory). The numerical objects that one has to deal with are arrays of size M"^^^"^ 
and M^^ . In practice, the calculations have been performed on conventional IBM RS/6000 workstations (43P-260 
and 3CT). The computation time for an iteration over a typical time span [0,40] depends strongly on the chosen 
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parameters. It ranges from several milliseconds for M = 3, K = 2 (program size 7 MB) over several seconds for 
M = 3, K = A (program size 8 MB) up to several hours for M = 5, K = A (program size 176 MB). The strongly 
limiting factor is the program size since the size of the arrays grows exponentially with K. E.g., the parameter 
combination M = A, K — 6 leads to too large arrays and cannot be treated by standard programming techniques. In 
practice, with the choice M = 6, K = 3 or M = 5, K = A, we already are at the upper limit of the QUAPI algorithm. 



We proceed in reporting our results for the specific example of the parametrically driven dissipative quantum 
harmonic oscillator. With the reduced density matrix (^7|) at hand, we can calculate the variances ( p^ ) within the 
QUAPI algorithm and compare them with the analytic predictions (0|9|). Most of the figures contain results for 
rather extreme parameter values, e.g. low temperature and large driving amplitude, in order to show that the QUAPI 
algorithm performs satisfactorily also in these limits. For more moderate choices of the parameters, the agreement 
(not shown) between numerical and analytic results is much better. 

Our main goal is to study the dependence of the variances (l^ on the QUAPI parameters M, K and At. For finite 
M and K, the deviation increases proportional to At^ due to the Trotter splitting in Eq. (|2^) with increasing At. 
For decreasing At, the Trotter error decreases but the error made by the memory truncation in Eq. ( p^ ) starts to 
dominate since more and more bath correlations are neglected. Thus, the overall error increases again. In between 
there exists an "optimal time step of least dependence" , where the quantities are least sensitive to variations of At. 
This represents the "principle of minimal sensitivity" for the optimal choice of the time step At for the QUAPI 
algorithm (see also [^). For M finite and K — > oo, the result would be independent of At for small At since the 
Trotter error would vanish and also the finite- memory error would not exist. 

The choice of M and K should be adapted to the chosen bath parameters. In the case of no driving, if the 
temperature is low, only few energy eigenstates are required, i.e. M may be chosen small. However, low temperature 
induces long-range bath correlations. Therefore, the memory length K has to be assumed large. The opposite holds 
true in the other limit of high temperature. In the case of driving, the number M of basis states is more important 
compared to the undriven case, since the variances oscillate strongly and higher lying energy states are excited. The 
memory length K has to be reduced instead if one is interested in the oscillation amplitudes. However, for the mean 
value of the variances, the total memory length K is again more important and should be maximized (see below). 



We shall choose two representative parameter sets for our considerations. Since the memory in Eq. ( |28| ) is truncated 
in the QUAPI algorithm according to Eq. (p9[), the crucial parameters are the temperature T and the damping 
strength 7. The relatively high temperature T — l.O and the small damping 7 = 0.1 form the first parameter 
set (High temperature - weak damping). For this choice, the numerical results are expected to agree well with the 
analytic results because large T suppresses the long-time memory contributions in Eq. ( p8[ ) and additionally, a small 
7 diminishes the influence of the bath correlations (|l^). Our second parameter set is given by T — 0.1,7 — 1.0 (Low 
temperature - strong damping). In this case, long-range bath correlations ( p^ ) play a major role and the truncation 
of them will induce an error which will be larger than in the case of a high temperature and weak damping. For 
intermediate parameter regimes, we find no qualitative differences. 

In all our calculations, we set to = and choose as the initial state the ground-state of the maximally curved (i.e 
— > + 2e) harmonic oscillator, i.e. ps{to = 0) = |0)(0|. The corresponding initial variances in Eq. ( |l7|) readily 
follow as cr^^ — l/{2y/u!Q + 2e), a^^ — and a-^p = lUq + 2e/2. Our standard choice for the cut-off function will 
be the exponential cut-off (^, if nothing else is stated. Furthermore, we always choose the diniensionless curvature 
loq = 1.0 in order to have a rather small separation of the energy levels in the undriven oscillator. This induces a high 
sensitivity on the number M of basis states since the higher lying states are then easily populated thermally or by 
driving induced transitions. The choice of a larger loq would be more in favour of the numerical algorithm. 



First, we consider the undriven case e = 0. Fig. |l| depicts the results for a high temperature T = 1.0 and small 
friction 7 = 0.1. Here and in the following, we use the dimensionless quantities which have been introduced in Eq. 
(|l^). Moreover, luq — 1.0 and ujc = 50.0. We find very good agreement with the analytic solution for the variances. 
The initial transient oscillations are reproduced and the asymptotic values for long times as well. The initial jump of 
(^xp{t) (of Eq. (pisi)) is still visible, while the jump of (Jppit) is proportional to 7^ and is not visible on this plot. 

To be able to study the dependence of the QUAPI algorithm on the parameters A/, K and At we consider the 
asymptotic values of the variances at long times. It is clear from Eq. ( p9| ) that (7xp{oo)—Q, so we focus in Fig. || on the 



IV. RESULTS 




A. High temperature — weak damping (no driving) 
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two non-trivial variances (Txxioo) and tTpp(oo). The qualitative dependence of both variances on the time step At is 
always similar: The deviation increases with increasing At due to the error proportional to At^ in the Trotter splitting 
in Eq. (p5|). For decreasing At, this error decreases and the "finite-K" -error takes over. The relevant At- value on 
which we focus in the following is the one for which the numerical result varies the least, i.e. the minima in the curves 
in Fig. ^ ("principle of minimal sensitivity" |T^]). 

The left column of Fig. || confirms that for a fixed memory length At ■ K, a smaller time step At induces a smaller 
Trotter-error whereas the finite-K-error remains roughly the same. While for a fixed M (left column in Fig. ||) QUAPI 
tends to underestimate the analytic result as K increases, a fixed K and growing AI (right column) leads to an opposite 
trend, suggesting that indeed the analytic result will be approached best when both M and K become large (at a 
plateau At- value tending towards zero). 

B. Low temperature — strong damping 

1. No driving 

Fig. |3| depicts the time-dependence of the variances in satisfactory agreement with the analytic result. The initial 
jumps of (Txp{t) and of (Jpp{t) are more pronounced in this strong damping case since the jumps are proportional to 7 
and 7^ (see Eq. (p^). The deviations in the transient behavior are due to the assumption of a strictly Ohmic classical 
dynamics (infinite cut-off ujc) in the analytic solution, see the discussion at the end of section They become more 
pronounced for low temperatures and strong friction because this assumption induces deviations in the short-time 
evolution of the variances on a time-scale The bath- induced long-range memory at this low temperature carries 
the deviations over the whole range of the transient dynamics. The fact that the memory length K is decisive for this 
low temperature is confirmed by the dashed-dotted line. 

The dependence of the asymptotic values axx{oo) and (Tpp{oo) on the QUAPI parameters is shown in Fig. ^. The 
number M of basis states is not so important, while the memory length K is decisive. Again, the analytic prediction 
is correctly approached when both M and K are increased. 

2. With driving 

Fig. H demonstrates for a small driving amplitude reasonable agreement with the analytics. The long-memory 
parameter set with K = 6 hits best the asymptotic mean value, but the oscillation amplitudes and frequencies are 
obtained best by the choice of a large M = 5. In comparison to the undriven case (e = 0) the time averaged variances 
are almost unchanged (Figs. |4|and ^) while the time-resolved behavior (Figs. ^ and ^) displays notable differences. 

Fig. depicts the time evolution for the relatively large driving amplitude. As expected, for strong driving, a large 
number M of basis states are required to describe the oscillations correctly. The averaged asymptotic values axx{oo) 
and app{oo) are plotted in Fig. |^. Since the strong driving mixes high energy eigenstates, the results are considerably 
more sensitive to the choice of M than for weak driving (Fig. ^ upper right panel). However, the same argumentation 
applies like in the undriven case (see Fig. ^). Considering the rather extreme parameters (small level-spacing, strong 
driving, low temperature, strong damping) the agreement with the analytic results is still satisfactory. 

C. Diverging dynamics and dependence on the cut-off ujc 

Fig. ^ shows Uxxit) for parameters belonging to an instability region of the Mathieu oscillator ( p^ ) p8[ , i.e. the 
variances for the driven quantum harmonic oscillator diverge for long times. Since the QUAPI algorithm is restricted 
to a (finite) M-dimensional Hilbcrt subspace it cannot reproduce such an asymptotic divergence. 

The last issue we address is the dependence of the dynamics on the cut-off parameter lOc and on the explicit shape 
of the cut-off function (^,^. First, we keep an exponential cut-off but choose a smaller cut-off frequency ujc- It 
is well known that for the (undriven) quantum harmonic oscillator cTpp (00) diverges with uJc, while axx{oo) is 
asymptotically independent of cuc- In Fig. |lO|, we choose the "worst" case (i.e. low temperature and strong damping) 
without driving and decrease the cut-off to uJc — 10.0. Compared to Fig. 0, the value of (Txx{oo) is indeed practically 
unchanged while cTpp (00) has notably decreased. 

Fig. |ll] shows results for a step-like cut-off (||). First, we observe that mainly the short-time behavior of the 
relaxation process is affected. Clearly, QUAPI with its restriction to only a few energy eigenstates cannot reproduce 
the transient high-frequency oscillations of (Jpp{t). Second, we note that a step- like cut-off affects the decay of the 
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response function L{t) from Eq. (|T^) for t oo. The real/imaginary part of L{t) decays qualitatively like an 
algebraically damped cos/sin- function. While this might suggest a strong dependence of the QUAPI results on the 
memor y le ngth, we actually find a rather weak dependence since the agreement between numeric and analytic results 
in Fig. is not considerably worse than in Fig. |. This means that the memory truncation in Eq. ( p9|) is in fact not 
very sensitive to the choice of the cut-off function /c(w, uJc) as long as one is not interested in the detailed short-time 
behavior. 



V. CONCLUSIONS 



We have studied the dependence of the QUAPI algorithm on its three numerical parameters, namely the time-step 
At, the number M of basis states, and the memory length K. As a test system, we have used the analytically solvable 
dissipative quantum harmonic oscillator and its parametrically driven generalization. The comparison shows a decent 
agreement of the approximative numerical result with the analytic solution, even in the case with driving. This means 
that a spatially continuous system can be described reasonably well by taking only a few basis states and a finite 
memory length into account. For low temperatures and weak-to- moderate driving, the number M of basis states has 
to be chosen small and the memory length K large, while in the opposite regime of high temperature, M has to be 
large but K may be chosen small. In both cases, satisfactorily large M and K values are still numerically feasible. 
For strong driving, the deviations increase but the QUAPI results are still in qualitative agreement with the analytic 
predictions. 

Our findings demonstrate the reliability of the QUAPI algorithm even in driven, spatially continuous systems 
and not only in finite, discrete dissipative quantum systems such as the spin-Boson-system. Therefore, the QUAPI 
algorithm may become a standard procedure for simulating open quantum systems in the presence of a novel class 
of time-dependent, not necessarily periodic driving fields. This technique is especially interesting for the study of 
decoherence in interacting two-level-systems processing quantum bits. There, the quantum gate operation prescribes 
the time-dependence of the external control fields which may exhibit a complex non-periodic time-dependence. 
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FIG. 1. Time-dependence of the variances axx(t), cjxpit), and o-pp(t) for the undriven dissipative quantum harmonic oscillator 
(e — 0,uio = 1.0) with bath parameters T — 1.0,7 — 0.1 and an exponential cut-off Eq. (^ with ujc — 50.0. In all the figures, 
we have used dimensionless quantities according to Eq. (|l^). The solid lines depict the analytic results (|l^) while the dashed 
lines represent the numerical solution obtained by the QUAPI algorithm with AI — 5, K = 4, At — 0.2. The asterisks mark the 
initial variances cr°^ = a^p — 0.5 and a^p = 0. 
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FIG. 2. Asymptotic values of the position (upper row) and momentum (lower row) variances Oxx{oq) and Opp{oo), respec- 
tively, as a function of the time-step At and different combinations of the two QUAPI parameters M (number of basis states) 
and K (number of memory time steps) for the undriven dissipative quantum harmonic oscillator (e = Q,ujo = 1.0) and heat 
bath parameters T — 1.0, 7 = 0.1 and uJc = 50.0 with exponential cut-off Eq. For the left column figures, the number M 
of basis states is fixed to M = 5 and the memory length K is varied, while for the right column figures, K is fixed to K = 'i 
and M is varied. Interconnected symbols: solutions obtained by QUAPI. Horizontal solid line: analytic result wjyi- 
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FIG. 3. Same as Fig. |l|, but for the parameters T — 0.1,7 — 1.0. Here, the QUAPI parameters are M = 2i,K 
(dashed line) and M = 5, iC = 4, At = 0.2 (dashed-dotted line). 
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FIG. 4. Same as Fig. but for the bath parameters T = 0.1, 7 = 1.0. 
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FIG. 5. Time-dependence of the variances i7xx{t),i7xp{t) and o-pp{t) for the parametrically driven dissipative quantum har- 
monic oscillator with luq = 1.0 and a small driving amplitude e = 0.1. The bath parameters are T = 0.1,7 = 1.0 and ujc — 50.0 
(exponential cut-off Eq. (||)). The QUAPI parameters are M = 3, A" = 6, At = 0.2 (dashed line) and M ^ 5, K = 4, At ^ 0.25 
(dashed-dotted line). The asterisks mark the initial variances crS^ = 0.45, a°p = 0.55 and a°p = 0. 

tt)„=1 .0, £=0.1, T=0.1, 7=1.0,(0^=50.0 





G — oM=3 




H hM=4 


,A'' 


a---aM=5 




A A'' 




a - ^r^J^ 


K=4 





0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0 
At At 



FIG. 6. Time-averaged asymptotic values of the position (upper row) and momentum (lower row) variances Wxx{oo) and 
CTpp(oo), respectively, versus time-step At for different combinations of the QUAPI parameters M and K, small driving ampli- 
tude e = 0.1 and bath parameters T = 0.1,7 = 1.0, loc = 50.0 (exponential cut-off Eq. (^). For the left column figures, the 
number M of basis states is fixed to AI — 5 and the memory length K is varied, while for the right column figures, K is fixed 
to = 3 and M is varied. The oscillator frequency is uiq = 1-0. 
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FIG. 7. Same as Fig. but for the strongly driven case e = 0.5. Here, the QUAPI parameters are M = i, K = 4, At = 0.25 
(dashed line) and M — 5, K — 4, At — 0.25 (dashed-dotted hne). The asterisks mark the initial variances (jL = 0.35, (jOp = 0.71 
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FIG. 8. Same as Fig. but for the strongly driven case e = 0.5. 
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FIG. 9. Time-dependence of the position variance axx{t) for a parameter set where the classical dynamics is unstable, i.e. 
e = 0.5,7 = O-l- The temperature is T = 1.0 and = 50.0 (exponential cut-off Eq. (^)). The QUAPI parameters are 
M — 5, K = i, At — 0.25. The asterisk marks the initial variance cr^^ — 0.35. The oscillator frequency is ujq = 1-0. 
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FIG. 10. Time-dependence of the variances axx{t),axp{t) and crpp(t) for a small cut-off frequency ujc = 10.0 (exponential 
cut-off Eq. (|)), LJo = 1.0, e = 0,7 = 1.0 and T = 0.1. Here, the QUAPI parameters are M = 5,K = 4,At = 0.30. The 



asterisks mark the initial variances a'. 
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FIG. 11. Same as Fig. but for a step-like cut-off Eq. (g) with lJc = 50.0. Parameters are ljq — 1.0, e = 0,T — 0.1 and 
7 = 1.0. The QUAPI parameters are AI = 5, K — 4, At = 0.25. The asterisks mark the initial variances a'^^ — — 0.5 and 



a^p = 0. 
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